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ABSTRACT 

A study was conducted of the digital computer simulation 
of the equations of motion of the Surface Effect Ship with 
six degrees of freedom as developed by the Oceanics Corpora- 
tion. The version of the program under consideration was a 
simulation of the Bell Aerospace Systems 100 -B Captured Air 
Bubble Craft (CAB) as adapted to the IBM-36O/67 digital 
computer and operating in irregular seas. The objective of 
the study was to optimize the simulation by reducing compu- 
tation time required to obtain solutions to the forces and 
moments acting on the craft while maintaining a reasonable 
degree of accuracy of the output variables. CPU time 
savings achieved by use of the FORTRAN H compiler are 
documented. The dependence of CPU time and output accuracy 
on the tolerance levels established for the integration 
algorithm is discussed. CPU time savings versus output 
accuracy as the tolerance levels are increased is presented. 
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TABLE OP SYMBOLS AND ABBREVIATIONS 



INTGRL 


= Integration Subroutine 


RHS 


= Subroutine for calculating the right hand side 
of the differential equations which describe 
the SES 


F 


= Force (Subscripted to indicate direction) 


g 


= Gravitational constant 


Lx 


= Mass moment of inertia about x-axis 


Lz 


= Mass product of inertia in x-z plane 


Ly 


= Mass moment of inertia about y-axis 


I 

zz 


= Mass moment of inertia about z-axis 


K 


= Roll moment 


m 


= Mass 


m b 


= Mass of air bubble 


P 


= roll angular velocity 


P b 


= pressure in the bubble 


Q 


= Pitch angular velocity 


Q in 


= Volumetric flow rate into plenum 


^out 


= Flow rate due to leakage 


R 


= Yaw angular velocity 


U 


= longitudinal velocity 


V 


= lateral velocity 


V/ 


= vertical velocity 


X 


= horizontal distance in direction of forward motion 


Y 


= horizontal distance positive to starboard 


Z 


= vertical distance positive downward 
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pitch angle, positive bow up 
standard atmospheric density 
roll angle, positive port side up 

yaw angle, positive turn to starboard 

d( ) 

dt 
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I. INTRODUCTION 



A. BACKGROUND 

In recent years the U. S. Navy has developed an interest 
in the Surface Effect Ship (SES) concept for combatant and 
support type craft. As a result of this interest, Oceanics 
Incorporated was commissioned by the U. S. Navy to develop 
a digital computer simulation of the six degrees of freedom 
equations for SES craft that would yield time domain outputs 
of loads and motions. This initial report was delivered to 
the Navy Surface Effects Ships Project Office (SESPO) and 
dated August, 1971 [Ref. 1], 

The simulation program was designed to achieve a suffi- 
ciently accurate and stable solution to the loads and motions 
equations for the modeling of the craft. The desire for real- 
time solutions was acknowledged, as they would be valuable 
for manned control simulation studies and greatly reduce 
costly computer time during case studies [Ref. 1], 

The program was developed using FORTRAN IV computer 
language and was prepared in modular form for ease in adapta- 
tion to various SES configurations. The program delivered 
to the Naval Postgraduate School in October, 1972, was a 
simulation of the Bell Aerospace Systems 100-B Captured Air 
Bubble Craft (CAB), a 100 ton craft. 

Results of simulated runs conducted by Oceanics, 
utilizing a Control Data Corporation 6600 digital computer 
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(using the scope 3-3 operating system) indicated that 
problem time to computer time ratios, as low as 3:1 for 
calm water runs and as large as 1:4 for large amplitude, 
oblique, irregular sea cases at high speed, could be 
expected [Ref. 2]. 

The initial studies at the Naval Postgraduate School, 
using the IBM-360/67 computer with the standard Fortran G 
compiler to solve the 100-B loads and motion program for 
sea state conditions, resulted in problem time to computer 
time ratios of 1:40 to 1:70, depending on the type run 
simulated. As the sea state condition increased, the 
computer time increased to the point where for certain type 
runs-, 2 to 3 hours of computer time were needed for problem 
solution [Ref, 3]« For example, broaching studies simulating 
65 knots with sea state 3 (irregular waves) at 150° encounter 
angle and a 30 second run time, required 121 minutes of CPU 
time. It should be noted however, that in the broaching 
study, the output requirements included a listing of sidewall 
gap each time subroutine RHS is entered thus adding a signif- 
icant amount of CPU time. 

B. OBJECTIVES 

The high ratio of problem time to computer time for sea 
state operations demanded that an in depth study be made of 
the Oceanics SES simulation program as applied to the IBM- 
360/67 located at the Naval Postgraduate School's W.R. Church 
Computer Facility with the following objectives in mind: 
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1. ) Careful analysis of each subroutine to determine 
those areas where the bulk of the computer time Is being 
spent . 

2. ) Determine methods whereby the computational effi- 
ciency of the program, can be maximized with minimum loss 
In accuracy. 

This thesis examines the SES 100-B digital computer 
simulation program on the IBM- 36 O /67 computer with the 
above mentioned objectives In mind and makes recommendations 
for substantially reducing computer time while maintaining 
a reasonable degree of accuracy with regard to the loads 
and motion simulation. 



II. INTEGRATION METHOD 



A. INTRODUCTION 

The SES simulation program is based on a mathematical 
description of the craft in six degrees of freedom as 
reported by Kaplan, Bentson, and Sargent [Ref. 1]. The 
resulting equations, in a form suitable for numerical 
integration on a digital computer, are given by 
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X = u 
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y = v 



(ID 



ip = R 



( 12 ) 



^a^in ^out^ 



(13) 



where X, Y, and Z, are the 3-dimensional spatial coordinates; 
P, Q, and R, are the angular velocities about the X, Y, and 
Z, axis respectively; I is the mass moment of inertia, and 
P is the force. is the mass of bubble air in the plenum; 

Qin Is the volumetric quantity of air flow rate into the 
bubble; Q qu ^ is the flow rate due to leakage, and p & is the 
standard atmospheric reference density. 

The program is of modular design utilizing a parti- 
tioning technique on the craft to calculate the various 
forces and moments exerted on the hull as it moves through 
the water. The forces and moments acting on the craft are 
computed in various subroutines using initial conditions 
supplied by subroutine INCON. Subroutine RHS consists of 
a simple summing algorithm to compute the total of the forces 
and moments in six degrees of freedom. 

Subroutine INTGRL calls subroutine RHS for the calcula- 
tion of the right hand side of equations defined by the 
Runge-Kutta-Merson (RKM) algorithm for each required step 
size. The elements of each k vector in the RKM algorithm 
are the ten variables defined by equations (1) through (9) 
and equation (13). 



16 



Since the variables X, Y, and ip defined in equations 
(10), (11), and (12), are expected to have slow variation, 
it was decided that a simple means of integration using 
Simpson's Rule could be applied to these equations. These 
variables are calculated in the main program [Refs. 1 and 2]. 

A complete description of the computer program including 
its organization, flow charts, input and output description, 
program listing, and a user's manual for the operation of 
the program are available in Refs. 2 and *J. 

B. SUBROUTINE INTGRL 

The numerical integration technique used in subroutine 
INTGRL is the Runge-Kutta-Merson (RKM) method for solution 
of a system of first order differential equations of the 
form 



y = f(t,y) 



(1*0 



the recursion formula is 




( 15 ) 



where 




( 16 ) 




( 17 ) 
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( 18 ) 
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( 19 ) 



k 5 3 f + k + 2 ~ 2 ^ + ^ 



( 20 ) 



and h is the stepsize. The truncation error in ( 15 ) is 
given by 



Initial conditions are read in through subroutine INCON 
to establish the starting point for each of the ten variables 
plus a trial stepsize (h) and a maximum error tolerance level 
for each integrator. Computations are then carried out and 
the truncation error given by equation (21) is computed. If 
any element in the error vector e is greater than the corre- 
sponding maximum tolerance level then the trial time stepwise 
is halved and the computations are repeated. If the error is 
less than the tolerance level, but greater than 1/16 this 
level, the computations proceed in accordance with the RKM 
algorithm. If all elements in e are less than or equal to 
1/16 the maximum tolerance level, the stepsize is doubled 
in the next integration step. This procedure automatically 
adjusts the time step in the computation to reflect the 
nature of the perturbations being experienced by the craft. 




( 21 ) 
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As a result, for slowly-varying phenomena (e.g., calm 
water, slow speed simulation), the computations are carried 
out quickly since larger time steps are taken, while high 
frequency effects (e.g., heavy seas, high speed simulations), 
will cause an increase in computation time. 
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III. PROCEDURE 



A. INTRODUCTION 

The obvious starting point for an efficiency study of 
the SES simulation program was to gain a thorough knowledge 
of the force and moment algorithm used. The object was to 
study the interconnections and interactions between the 
various subroutines with a goal of arriving at a more effi- 
cient (i.e. less time consuming) solution by reprogramming 
and/or better utilization of the IBM -360 software. 

A major reprogramming effort was considered to be beyond 
the scope of this thesis. However, investigation of subrou- 
tines which were known to consume the bulk of the computation 
time was conducted. In conjunction with this investigation 
various timer functions built into the IBM- 36 O were utilized 
to pin point time consuming procedures. 

B. IBM-360 COMPILER 

The SES simulation program was originally adapted to the 
Naval Postgraduate School IBM-360 computer utilizing the 
standard FORTRAN G compiler. This compiler features a DEBUG 
facility and core optimizations which lends itself to the 
smaller job, teaching environment prevalent at the school 
[Ref. 5]. Among other compilers, the IBM-360 computer system 
has available FORTRAN H, an IBM processor which produces a 
highly optimized code which is highly desirable for large 
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production programs such as the SES simulation program 
[Ref. 6 ], 

Compilation of the program In FORTRAN H resulted In an 
average CPU time savings of thirty-seven percent over the 
FORTRAN G compiler while yielding exactly the same output 
variable values. For example, a simulation of a thirty 
second run in a state three, head sea compiled in FORTRAN G 
required 80.5 minutes of CPU time. The same simulated run 
time compiled in FORTRAN H required only 50.7 minutes of 
CPU time. A number of runs, under varying initial conditions, 
were made to verify the CPU time savings using the FORTRAN' 

H compiler. Unless otherwise indicated all further runs 
discussed in this thesis were made using the FORTRAN H 
compiler . 

C. TIMER SUBROUTINES 

The IBM- 36 O /67 at the Naval Postgraduate School’s W.R. 
Church Computer Facility has several built in timer subrou- 
tines which were helpful in determining those areas in the 
program in which the bulk of the CPU time was being used. 

Two such routines were used in this study. 

1 . Subroutine IONUM 

Subroutine IONUM is a program which monitors the 
number of times an input or output device is accessed during 
the run of a computer program. In adapting the Oceanics 
simulation program to the IBM-360, extensive use of discs 
as temporary storage devices was made. By monitoring the 
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number of input/output operations with IONUM a more effi- 
cient data blocksize for the discs could be determined, 
thus reducing the number of accesses required. 

2 . Subroutine PROGLOOK 

Several areas of the SES simulation program were 
suspected to be inefficient and/or requiring excessive time 
during the program solution. In order to pin point these 
areas of the program, subroutine PROGLOOK was used. 

PROGLOOK, when linked to the program under consid- 
eration, monitors the sequential flow of the program instruc- 
tions. The routine "strobes" the instruction flow at a rate 
of 50 samples per second or approximately once every 4000 
instructions. In this manner a time history of the program 
functions was obtained. The output of subroutine PROGLOOK, 
which was of interest to this study, is in the form of a 
histogram which gives a plot of percentage of time the program 
spent performing various functions. Correlation of the loca- 
tion of these functions with a computer map listing of the 
SES simulation program served to point out time consuming 
areas within the program. 

D. INTEGRATOR TOLERANCES 

A prime suspect area of the simulation program for excess 
time consumption was the Runge-Kutta-Merson (RKM) algorithm 
used to solve 10 of the force and moment differential equa- 
tions. If during a given time increment, any one of the 10 
error function values is found to be outside the given 
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tolerance, the stepsize Is halved and all variables are 
computed using the new stepsize. This process is continued 
until all 10 error functions are within their respective 
tolerances and then a new time increment is started. 

The above stated procedure does not constitute an inord- 
inate amount of computer time under calm water simulations 
as the variable values tend to vary slowly and the RKM 
algorithm converges rapidly to the established tolerance 
values. However, craft perturbations caused by such condi- 
tions as simulations of rapid turns, irregular waves to 
simulate sea state, and loss of a propulsion unit, can cause 
some of the variables to vary at a high rate of change, thus 
requiring many iterations of the algorithm for convergence 
to the established tolerance levels and resulting in large 
CPU times. 

Reference 3 established the standard maximum tolerance 
levels for each integrator according to procedures set forth 
in Ref. 2. These tolerance levels were considered to be 
the basis of accuracy desired for the SES' simulation program 
used at the Naval Postgraduate School and are listed in 
Table I. The input and output variables for each of the 
13 integrators are shown in Fig. 1. 
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Table I 



Tolerance Values For Integrators 



Integrator Number 
1 
2 
3 
H 

5 

6 

7 

8 • 

9 

10 



Output Variables 
U 
V 
W 
P 
Q 
R 

4 > 

6 

Z 

"b 



Tolerance Values 
.0001 
.0002 
.000001 
.0000001 
.0001 
.000001 
.000000001 
.000000001 
.0001 
.0001 



1. Steady State Conditions 

A number of runs were conducted initially, for each 
simulation study made, in order to establish steady state 
conditions for a given set of initial conditions. Reasonably 
accurate steady state conditions are essential to prevent 
initial transients in the simulation program which can 
cause the program to abort, or worse yet, to give erroneous 
results. The key variables to initialize for a given input 
are thrust, draft, plenum pressure, and to a lesser degree 
pitch angle. 



2M 



, 




Q 




8 







RKM 





RKM Runge-Kutta-Merson Method 
SR Simpson's Rule 

Figure 1. Inputs and Outputs of Integrators 
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Two run conditions were selected to fulfill the 
objectives of this thesis. (1) Sea state 3, 40 knot speed, 
head seas, and (2) Sea state 3, 60 knot speed with turn. 
Table II gives the initial conditions for these two run 
conditions. 

Table II 

Initial Conditions 
(Sea State 3) 



Speed 
(kts ) 


Draft 

(In) 


Plenum 

Pressure 

(PSF) 


Thrust 
(one Engine) 
(Lbs) 


Pitch 

Angle 

(Deg) 


40.0 


30.0 


78.0 


7000.0 


O 

• 

VJ! 

o 


60.0 


28.0 


80.0 


10800.0 


0.45 


2. 


Integrators 


3 and 10 







The greatest effect of the sea state condition on 
the SES program computation occurs on integrator 3 (heave 
accelerations) and integrator 10 (time rate of change of 
the bubble air mass). The majority of the forces and 
moments acting on the craft when encountering sea states 
is manifested in changes in Z directed accelerations and 
plenum volume (i.e., bubble mass). 

As a first approach to reducing the CPU time for 
sea state simulations, the changing of the implementation 
of the RKM algorithm was considered. It was felt that if 
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Integrator 3 and 10 could be Isolated and their solutions 
made to converge before the other eight integrations were 
attempted, then a considerable computation time savings would 
be gained due to the reduced number of iterations needed 
through subroutines INTGRL and RHS . However,. due to the 
functional dependence of heave acceleration and bubble mass 
on some of the variables, it was found that Integrators 3 
and 10 could not be made to converge independently of the 
other eight. 

Since both Integrator 3 and 10 were extremely sensi- 
tive to craft perturbations the decision was made to loosen 
the tolerances on these two integrators and compare the 
results with those obtained using the standard tolerances 
(See Table I) . It was felt that if the percentage change 
in variable values was not too great (less than 10 percent), 
acceptable results could be obtained with a considerable 
CPU time savings. 

The two integrator tolerance levels were increased, 
in tandom, by a factor of ten in steps of one for a total 
of ten runs under a given set of initial conditions. Each 
computer run simulated a 30 second craft run and parameter 
values were printed out every 0.05 seconds for a sample 
rate of 20 per second per variable. In addition, the root 
mean square (RMS) value and the average value of the 
parameters was calculated over the entire run interval for 
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comparison with values obtained from runs using the standard 
tolerances. Using these results an error analysis could 
be conducted. 

3 . Integrators 1,2 and 4-9 

Analysis of the CPU time versus output accuracy 
as observed in the above runs indicated an optimum tolerance 
of 0.000007 for integrator 3 and 0.0007 for integrator 10. 

It was then decided to investigate the sensitivity of the 
other eight integrators. The method used for this investi- 
gation was to hold the above tolerances on integrators 3 
and 10 while increasing one of the other integrator toler- 
ances by a factor of ten (the other seven tolerances remain- 
ing at standard values). Each integrator was tested in 
succession and the results compared to standard runs. 
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IV. DISCUSSION AND EVALUATION OP RESULTS 



A. FORTRAN H COMPILER 

The efficiency of the FORTRAN H compiler can best be 
illustrated by the following example. Reference 3 cited 
the disparity in CPU times between simulation runs using 
control statement 00201 and 00202. These runs were made 
using the FORTRAN G compiler. Control statement 00201 inputs 
directly to subroutine INCON the overall weight of the craft 
plus the location of the center of gravity and the mass 
moment of inertia about the X, Y, Z, and XZ-axis. Control 
statement 00202 inputs distributed weights and moments and 
subroutine INCON must then compute the craft weight, CG 
and moments. It would appear that the additional computa- 
tions involved in using control statement 00202 would require 
more CPU time than statement 00201. 

Results of Ref. 3 indicated that the opposite result . 
was obtained and that CPU time decreased by approximately 
10 percent when control statement 00202 was used as input. 

The reason for this result has not been established and 
should be looked into by future investigators. Nonetheless, 
results of runs compiled in FORTRAN H reduced the CPU time 
difference cited in Ref. 3 to less than 1 percent which 
illustrates the optimization capabilities of this compiler. 



29 



B. INTEGRATOR TOLERANCES 



The initial runs were made to verify the sensitivity 
of the Runge-Kutta-Merson algorithm to integrator tolerances. 
The runs simulated 40 knots speed, state 3 seas, head waves, 
30 second run, and confirmed that integrators' 3 and 10 
were the most sensitive to craft perturbations. The results 
of these runs were as follows. 

Table III 

CPU Time Comparison 







(Sea State 3) 




Integrator 


Tolerance 


CPU Time 


PERCENTAGE De 


No. 3 


No. 10 


(Minutes ) 


In CPU Time 


0.000001* 


0.0001* 


53.55 


Standard 


0.000002 


0.0002 


38.74 


27.66 


0.000003 


0.0003 


33.39 


37.65 


0.000004 


0.0004 


31.38 


41.40 


0.000005 


0.0005 


29.20 


45.47 


0.000006 


0.0006 


28.90 


45.73 


0.000007 


0.0007 


26.99 


49.56 


0.000008 


0.0008 


26.19 


51.09 


0.000009 


0.0009 


25.03 


53.26 


0.000010 


0.0010 


24.50 


54.25 



* Standard for comparison 
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1 . Integrator 3 and 10 



The tolerances on integrators 3 and 10 were loosened 
during the above stated runs while the other eight integra- 
tor tolerances were maintained at their standard values. 

Table III shows the results of these runs with regard to 
CPU time and Figure 2 is a graphical representation of the 
data. 

As the integrator tolerances were varied the RMS 
and average values of the output variables were compared 
with those obtained using standard integrator tolerances. 

The percentage deviations from standard values were computed 
in order to have a measure of output variable accuracy as 
a function of CPU time savings. Table IV tabulates these 
results for maximum tolerance error levels. 

As was expected, the highest percentage deviations 
occured with heave acceleration (integrator 3) and the 
variables concerned with air flow (integrator 10). All 
other variables had insignificant deviations (less than 1 %) 
from standard values. In general, the percentage deviation 
of the output variables increased almost linearly as the 
integrator tolerances were loosened, with the maximum 
deviation occurring for maximum tolerance level. 

Increasing the Integrator tolerance levels above a 
factor of seven did not greatly decrease the CPU time required 
for problem solution but did tend to decrease the accuracy 
of the affected output variables. For example, the initial 
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FIGURE 2 
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U DEVIATION FROM STANDARD VALUE (PERCENTAGE) 



Table IV 



Output Variable Maximum Deviation 
Prom Standard Values 

Output Variable Maximum Deviations from 

Standard Value (Percentage) 





RMS 


Average 


Draft 


0.04 


0.07 


Pitch Angle 


0.82 


0.00 


Plenum Pressure 


0.44 


0.29 


CG Heave Acceleration 


2.50 


0.00 


Surge Speed 


0.03 


0.02 


X-Displa cement 


0.04 


0.04 


Pan Power 


1.39 


3.33 


Air Flow Into Plenum 


0.05 


0.95 


Air Flow Out of Plenum 


0.76 


1.22 


Bubble Drag 


0. 6l 


0.27 


Pitch Rate 


0.00 


0.00 



tolerance increase of a factor of seven gives a CPU time 
savings of 26.56 minutes, while increasing the tolerance 
level to a factor of ten achieves only another 2.49 minutes 
savings in CPU time. At the same time the percentage 
deviation in the value of heave acceleration increases from 
1. 8 % to 2 . 5 %. 
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2 . Integrators 1,2 and 4-9 



With integrator 3 and 10 tolerances at 0.000007 and 
0.0007 respectively, each of the other integrator tolerances 
were increased in turn by a factor of ten. Simulated run 
conditions were as previously described and CPU times were 
compared. No appreciable CPU time savings was observed 
during any of these runs. The average CPU time was 27-59 
minutes for these runs. The percentage deviation in the 
values of the output variable was negligibly small, indi- 
cating that these variables were much less susceptible to 
the moments and forces caused by the simulated sea state 
conditions . 

C. DIFFERENT TYPE RUN CONDITIONS 

Various run conditions were simulated in order to study 
the effect of loosening the tolerances on integrator 3 and 
10 while maintaining standard tolerances on the other eight 
integrators. Each computer run simulated a 30 second craft 
run. 

Table V shows the results of these runs with respect to 
CPU time savings. 

Comparison of the results shown in Table III and Table 
V indicate that the percentage of CPU time savings is about 
the same regardless of the type run simulated. Therefore 
increasing the tolerance levels on integrators 3 and 10 
by a given amount will result in a predictable amount of 
computer CPU time savings. 



Table V 



CPU Time Comparison For 

Different Type Run Conditions 
(Sea State 3) 



Speed 


Integrator Tolerance 


Wave Encounter 


CPU 


Percentage 










Time 


Decrease in 


(knots) 


No. 3 


No. 10 


Angle 


(min. ) 


CPU tire 


40 


0.000001* 


0.0001* 


150° 


52.09 


Standard 




0. 000007 


0.0007 




24.73 


52.52 




0.000010 


0.0010 




23. 47 


54.94 


60 


0. 000001* 


0.0001* 


Partial Turn 


68.04 


Standard 




0. 000007 


0.0007 




34.34 


49.24 




0.000010 


0.0010 




31.93 


53.07 


60 


0.000001* 


0.0001* 


180° Turn 


' 69.97 


Standard 




0. 000007 


0.0007 




36.36 


48.03 




0.000010 


0.0010 




32.46 


53.61 



* Standard for comparison 

In comparing results of deviations of output variables 
under different run conditions, no change was observed 
between the 40 knot, head sea simulation and the 40 knot 
150° wave encounter simulation outputs. However, under the 
more severe conditions of 60 knots with turns applied some 
of the output variable values did have an increased deviation 
from the standard values. Table VI depicts these results 
for the variables concerned. 
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Table VI 



Output Variable Maximum Deviation 

From Standard Values 
(Sea State 3> 60 knots) 



Output Variable 



Plenum Pressure 



CG Heave 
Acceleration 



Fan Power 



Air Flow Out 
of Plenum 



Turn Condition 

Partial 

l80° 

Partial 

180 ° 

Partial 

l80° 

Partial 



Maximum Deviation From 
Standard Value (Percentage ) 



RMS 


Average 


1.55 


0.29 


1.30 


0.05 


4.76 


1.00 



5.13 


0.00 


0.90 


4.58 


1.41 


4 . 6l 


0.17 


1.69 



D. SUBROUTINE PR0GL00K 

In order to determine other areas of the SES simulation 
program where large amounts of CPU time was being used. 
Subroutine PR0GL00K was utilized. Integrator 3 and 10 
tolerances were 0.000007 and 0.0007 respectively. A 40 
knot, state 3 sea, head w aves, 30 second run was simulated. 
Figure 3 is a histogram of the results of the run. 

It was found that Hl% of the computing time was spent 
in subroutine V/ AVES. The bulk of this time was consumed 
in the wave table. A surprisingly large amount of time 
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OF CPU TIME SPENT IN VARIOUS 
SUBROUTINES AND FUNCTIONS 
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SUBROUTINE OR FUNCTION 
FIGURE 3 



(22%) was spent computing trigonometric function values. 

The bulk of the time spent in a given subroutine was gen- 
erally found to be involved in repetitious DO LOOPs. 

The area labeled "Other" included various FORTRAN built 
in functions and the remaining subroutines and functions 
of the SES program which were not identified on the graph. 

Each of these routines contributed less than 1% of computation 
time. 



38 



V. CONCLUSION AND RECOMMENDATIONS 



A. INTEGRATOR TOLERANCES 

The computer run time dependence of the SES simulation 
program on the tolerance levels chosen for the Runge-Kutta- 
Merson algorithm is clearly demonstrated by the results 
shown in this thesis. Further, CPU time is directly 
related to the magnitude of the tolerance levels of inte- 
grators 3 and 10 and is relatively independent of the 
tolerance levels established for the other eight integrators. 

Output variable accuracy as would be expected is also 
related to the tolerance levels of the integrators, however 
the percentage deviations from standard values is small and 
the values obtained could be considered adequate for most 
purposes. In particular, initial studies under a given set 
of conditions would realize a considerable amount of CPU 
time saving by using increased tolerance levels, while 
maintaining sufficiently accurate solutions to determine 
the relative magnitude of the forces and moments acting on 
the craft. As more precise values are required, the tolerance 
levels could be progressively tightened with an accompanying 
predictable increase in CPU time requirements (a valuable 
planning tool). 

It should be noted that the individual programmer should 
determine his standard of accuracy for each output variable. 
Since the magnitude of the outputs vary greatly, a small 
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percentage deviation in one value may be acceptable for a 
given case study while being entirely too large for another 
output value under the same conditions. For example, the 
air flow rate into the plenum has a magnitude on the order 
of 5000 cubic feet per second with a 1 % deviation from 
standard which may be entirely acceptable, whereas the CG 
heave acceleration magnitude of 0.32 "g" units with 1 % 
deviation may not be acceptable. 

Due to time limitations, the simulated runs conducted 
in this thesis were confined to relatively slow speeds in 
the 40-60 knot range. The deviations observed for the 
critical output variables of heave acceleration, plenum 
pressure, and air flow in and out of the plenum, increased 
with an increase in speed and with the introduction of turns. 
Therefore it is recommended that further studies be made at 
higher speeds and more severe sea state conditions to deter- 
mine the validity of results obtained with increased 
tolerance levels. 

B. DIGITAL PROGRAMMING TECHNIQUES 

Time studies of the SES Loads and Motion program using 
subroutine PR0GL00K indicated that an appreciable amount of 
computer time was spent in the trigonometric functions. Most 
large computers use nested Chebyshev Polynomials or Taylor 
Series expansions as an efficient and accurate solution to 
the trigonometric functions [Ref. 6]. A faster method of 
solution, for example, would be a four place table lookup 



for the functions. A careful analysis of the loss of 
accuracy incurred by this method would have to be conducted. 
Another consideration would be the increased core require- 
ments due to the addition of a lengthy table. 

It was noted that considerable computation time is spent 
in the execution of DO LOOPs throughout the program. It is 
recommended that an in depth study be made with the goal of 
replacing the iterative DO LOOP process with straight forward 
calculations where possible. An excellent example of this 
type of programming is presented in Appendix A of this thesis. 
A straight forward programming technique resulted in the 
elimination not only of DO LOOPs but a subroutine as well. 

C. HYBRID COMPUTATIONS 

Plans are being formulated at the Naval Postgraduate 
School with regards to the feasibility of adapting the SES 
Loads and motion simulation programs to a hybrid computation 
technique. The development of linearized equations of 
motion presented in Ref. 7, and the analog computer simulation 
developed in Ref. 8, could serve as a starting point for 
conversion to hybrid computation. 

The Naval Postgraduate School's hybrid computation 
facilities include the XDS 9300 digital computer interfaced 
with a COMCOR 5000 analog computer and an Adage Graphics 
Terminal. Perhaps the most formidable problem presented by 
a hybrid conversion would be the core storage requirement 
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in the digital computer. The present storage capabilities 
could however be augmented by disc and/or tapes. 

The parameter scaling problems referred to in Ref. 8 
could be readily handled in a hybrid configuration. Other 
advantages to hybrid computation of the Loads and motion 
program include near real time solutions and the versatility 
of being able to change parameters during the run. The 
visual representation of the output parameters afforded 
by the Adage Graphic Terminal would be a powerful engineering 
tool. Hybrid computation would be particularly useful in 
vertical plane studies such as heave acceleration analysis 
and prediction j and for control systems studies. 



APPENDIX A 



An Example of Programming Changes to Improve 

Efficiency 



Reference 9 presented the following changes to the SES 
Loads and Motions Program to improve the efficiency of the 
program. The affected subroutines are RHS and INCON and in 
addition the subroutine DMINV is deleted. 

The following portion of subroutine INCON was deleted: 

212 DO 211 I = 1,6 

DO 211 N = 1,6 

211 A(I,N) = 0.0 

DO 213 N = 1,3 

213 A(N,N) = AM 
A(4,4) = AIXX 
A( 5 ,5 ) = AIYY 
A( 6 , 6 ) = AIZZ 
A(4 ,6) = -AIXZ 
A(6 , *0 = -AIXZ 

AIMAX = AMAXI(AM, AIXX, AIYY, AIZZ, ABS(AIXZ) ) 

DO 214 I = 1,6 

DO 214 J = 1,6 

214 A( I , J ) = A ( I , J ) /AIMAX 
CALL DMINV(A,G,D) 

DO 215 I = 1,6 

DO 215 J = 1,6 

215 A ( I , J ) = A( I , J) / AIMAX 
IP (D.NE.0.0) GO TO 10 
WRITE (6,216) 

STOP 
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Inserted in place of the above programming was the 
following: 

AMASSI = 1.0/AM 
D = 1.0/(AIXX*AIZZ-AIXZ*AIXZ) 

DIXX = AIXX*D 
DIXZ = AIXZ*D 
DIZZ = AIZZ*D 
AIYYI = 1.0/AIYY 
GO TO 10 

The INC ON parameters were linked to subroutine RHS by 
the following COMMON statement 

y 

COMMON/ATRIX/AMASSI,AIYYI,DIXX, DIXZ, DIZZ 

In subroutine RHS the six element matrix GF(I) was 

deleted and the following identifiers were substituted for 

the summation of forces: SUMX, SUMY, SUMZ, SUMK, SUMM, 

SUMN. In addition the following deletion was made. 

DO 1 I = 1,6 
VALUE (I) = 0.0 
DO 1 J = 1,6 

VALUE (I ) = VALUE (I)+A'(I,J)*GF(J) 

1 CONTINUE 

Substituted for the above DO LOOPs was the following 

VALUE (1) = SUMX* AMASSI 
VALUE (2) = SUMY*AMASSI-R*U 
VALUE (3) = SUMZ*AMASSI 
VALUE (1*) = SUMK*DIZZ+SUMN*DIXZ 
VALUE (5) » SUMM*AIYYI 
VALUE (6) = SUMN*DIXX+SUMK*DIXZ 

The indicated changes deleted nine DO LOOPs and one 
subroutine resulting in a much more efficient program both 
from a time and a core requirement standpoint. 



APPENDIX B 



Equivalent Expressions for Yaw Acceleration 



Reference 9 included the classic Euler equations of 
motion in six degrees of freedom and listed the assumptions 
under which certain terms of these equations could be 
disregarded as being insignificant. The question was 
raised as to whether the equation describing Yaw Accelera- 
tion (R) in classic Euler terms was equivalent to the one 
defined by subroutines DMINV and RHS. 

This appendix shows that, under certain assumptions, the 
two expresssions are equivalent. 

The classic Euler equations as presented in [Ref. 9] are 
as follows. 



P 

U = — - QW + RY 
m 



(i) 



P 

V = — ^ - RU + PW 
m 



( 2 ) 



W = — - PV + QU 
m 



(3) 



P = 



F, I 
k zz 



F^ 1 
n xz 



Q I 



zz 



II - 1 
XX zz xz 



II -I 2 II -I 2 

XX ZZ XZ XX zz xz 



, P < I *x- I yy +I M ,I *« 
[ ^ 



- R(I„_ - I„„ + — ) ] 



'zz yy I 



zz 



(4) 
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( 5 ) 



A _ F m . p ««« - W 
Q " TZ + 13 + 1 



yy 



yy 



yy 



R = 



n 



PQ(I - I ) (P-QR)I 



zz 



XX 

T 



yy 



xz 



( 6 ) 



zz 



zz 



The fundamental equations of motion for six degrees of 
freedom developed in Ref. 10 are as follows: 

F x = m[U + QW-RV-X G (Q 2 +R 2 ) +Y g (PQ-R) +Z q (PR+Q)] (7) 

F y = m[V + RU - PW - Y q (R 2 +P 2 ) + Z q (QR-P) + X g (QP+R)'] (8) 

F z = m[W + PV - QU - Z q (P 2 + Q 2 ) + X q (RP-Q) + Y q (RQ+P)] (9) 

F k = I XX P+ (I zz -F yy )QR + m[Y G (W+PV-QU) - Z q (V+RU-PW) 

+ X g Y q (PR-Q) -X g Z q (PQ+R) + Y g Z q (R 2 -Q 2 )] (10) 

F rn = I yy < ^ + (I xx" I zz )RP + m[ z G (U+ Q W-R v ) - X q (W+PV-QU) 

+ Y q Z g (QP-R) - Y q X g (QR+P) +X q Z g (P 2 -R 2 )] (11) 

F n = I ZZ R+ (I yy - I xx )PQ + mCX G (V+RU-PW) - Y q (U+QW-RV) 

+ Z g X g (RQ-P) -Z g Y q (RP+Q) + Y g X g (Q 2 -P 2 )] (12) 

where X G , Y Q , and Z Q are distances from the center of gravity 
to the origin along the x, y, and z axis respectively. 
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Because of symmetry of the craft, the significant product 



terms of inertia are the XqZ q terms. Under these conditions, 
equations ( 7 ) — ( 12 ) reduce to the following form. 



P = m[U + QW - RV] 

X 


(13) 


F = m[V + RU - PW] 


(14) 


F = m[W + PV - QU] 
z 


(15) 




(16) 


P m * I yy« +(I xx- I zz )RP + I xz (p2 - R2) 


(17) 


P n - ***** 


( 18 ) 



where 



Z xz = mX G Z G 



Solving equations (13)- ( 1 8 ) for the desired derivatives 
yields : 



U = — - QW + RV 
m 



(19) 



V = -Z- - RU + PW 
m 



( 20 ) 
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( 21 ) 



W = 
P = 
Q = 
R = 



P 

-£■ - PV + QU 
m 



xx 



- ( 



I zz _I yy 

I 

xx 



P I -I 

m , xx zz 

_r yy _ 

F n ^yy'^xx 

' ' J zz 



)QR + I XZ 

XX 


(PQ+R) 


(22) 


I 


)RP - XZ 

yy 

i 


(P 2 -R 2 ) 


(23) 


) PQ ~ 


(RQ-P) 


(24) 



zz 



Equations (19)— (21) are identical to equations (l)-(3). 
Upon substituting equation (24) into (22) and rearranging 
equation (23) it can be shown that equations (22)-(24) are 
equivalent to equations (4)-(6). 

If it is assumed that P, Q, and R are each much smaller 
than one, then product terms involving them may be ignored 
and equations ( 1 ) — ( 6 ) reduce to the following equations as 
stated in Ref. 9- 



• 

U = 


P 

X 

m 


• 


(25) 


• 

V = 


P y 

TT - RU 




(26) 


w = 


P 

z 

m 




(27) 


• 

p = 


Vzz 


p I 
n xz 


( 28 ) 


JT 


II - 1 2 

XX zz xz 


+ 2 
II - I 
xx zz xz 
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(29) 



Q = 



m 

yy 



r = 



n 



zz 



+ P 



xz 



zz 



(30) 



The equations of motion as decoded from subroutines 
DMINV and RHS are as follows. 



F 

U = — 
m 



(3D 



V = - RU 



(32) 



W = z 



F 
2 

m 



(33) 



P = 



Q = 



R = 



f,i 

k zz 



FI 
n xz 



-r ^ 

II -I II -I 

xx zz xz xx zz xz 



m 



yy 



FI 
n xx 



P ^ 1 

k xz 



? 2 

II -I II -I 

XX ZZ XZ XX zz xz 



(3*0 



(35) 



(36) 



Equations (25)— (29) are identical to equations (31)-(35) 
Substituting equation (28) into equation (30) yields: 



1)9 



• F n P k I zz 

R = =£- + ( — 5 - 

i zz II -I d 

XX zz xz 



FI I 

n xz N xz 

T' T 

II - 1 zz 

XX zz xz 



F F. I 

= n . k xz 



F I ‘ 
n xz 



I zz I I - I ? 

XX ZZ xz 



I I 2 - 1 I 2 
XX ZZ zz xz 



F. I 
k xz 



R = 



I I 


-x 2 


XX zz 


xz 


F. I 




k xz 




I I 


C\J 

H 

1 


xx zz 


XZ 


F. I 




k xz 





+ F ( J- + xz 



n I tt 2 tt 2 

zz I I - I I 

xx zz zz xz 



I I 2 -I I 2 +I I 2 
X -C / XX zz zz xz zz xz N 

+ ±* s 2 ) 

I 2 (I I - 1 2 ) 

ZZ v XX ZZ xz 



F I 
n xx 



? ? 

II -I II -I 

xx zz xz xx zz xz 



(37) 



Equation (37) is identical to equation (36) therefore, based 
on the given assumptions, the classical Euler equations 
given in Ref. 9 and as developed from Ref. 10 are mathe- 
matically equivalent to those defined in subroutines DMINV 
and RHS. However it is felt that the validity of disregarding 
product terms should be more thoroughly investigated. It 
is recommended that equations of the form (l)-( 6 ) be 
incorporated in future digital simulation studies to determine 
the degree of contribution of the cross product terms on the 
resultant output variables. 
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